home *** CD-ROM | disk | FTP | other *** search
/ NetNews Offline 2 / NetNews Offline Volume 2.iso / news / comp / lang / c-part1 / 9362 < prev    next >
Encoding:
Text File  |  1996-08-05  |  4.1 KB  |  109 lines
  1. Path: zac001p06.zone.ca!user
  2. From: michel_saucier@cmq.qc.ca (Michel Saucier)
  3. Newsgroups: comp.lang.c
  4. Subject: Re: Speed of random generators
  5. Date: Sat, 09 Mar 1996 12:56:41 -0400
  6. Organization: Bell Global Solutions
  7. Distribution: world
  8. Message-ID: <michel_saucier-0903961256410001@zac001p06.zone.ca>
  9. References: <TAKAHASI.96Mar1113825@poisson.ece.cmu.edu>
  10. NNTP-Posting-Host: zac001p06.zone.ca
  11. X-Newsreader: Yet Another NewsWatcher 2.2.0b4
  12.  
  13. In article <TAKAHASI.96Mar1113825@poisson.ece.cmu.edu>,
  14. takahasi@poisson.ece.cmu.edu (Eduardo S. C. Takahashi) wrote:
  15.  
  16. > I'm working in a simulation that requires intensive use of pseudo-random
  17. > sequences. My problem is that profiling the program I found out that the
  18. > generation of the sequences is responsible for about 50% of the total 
  19. > execution time. Since each run can take about 4 hours, it would be good
  20. > if I could save some time. 
  21. > 
  22. > Presently, I'm using the function erand48(), and I'm also making some 
  23. > tests with ramdom(). If anyone has a suggestion I would really appreciate!
  24.  
  25. Here is the fastest pseudo-random number generator this side of chaos.
  26. It compiles to about 20 machine instructions, all 1-clock type, except for
  27. [cached] memory references, which take 1 or 2 clocks, depending on scheduling.
  28. All in all, alea() takes 30 clock cycles or less to generate a good
  29. quality 32 bits
  30. pseudo-random number. Don't take my word for it. Test it.
  31.  
  32. As a comparison, in most congruential generators (the most common type of RNG)
  33. there is two MULT (25 clocks) and two DIV instructions, each one costing
  34. about 40 clock cycles,
  35. for a total of about 140 clocks by value generated.
  36.  
  37. /* alea(): an halfway decent random number generator */
  38. static unsigned long AInd=0;
  39. static unsigned long AMem[64]; /* 32 values are needed, extended to 64 for
  40. a speedup */
  41. unsigned long alea(void)       /* LFSR-based random generator, with a
  42. 2*32=64 longwords memory */
  43. { unsigned long *p=AMem+(AInd++%32);  /* 32 LFSR are maintained in
  44. parallel, bitwise */
  45.   return *p=*(p+32)=*p^*(p+1)^*(p+2)^*(p+4)^*(p+6)^*(p+31);
  46.  /* each value is kept twice in the array, to speed up access */
  47. }
  48.  
  49. long ralea(unsigned long range) /* return a pseudorandom value in 1..range */
  50. { return 1+alea()%range;}
  51.  
  52. static long LHseed=1; /* please kindly CHANGE ME */
  53. unsigned long lehmer(void)/* Lehmer algorithm, Schrage's method to avoid
  54. overflow */
  55. {long hi=LHseed/127773L,lo=LHseed%127773L;
  56.  long t=lo*16807L - hi*2836L;
  57.  return LHseed=(t>0)?t:t+2147483647L;
  58. }
  59.  
  60. void initAMem(void) /* initialize the AMem[0..31] array columnwise */
  61. { int i,j; unsigned long t; 
  62.   LHseed=(unsigned long)clock(); 
  63.   for(i=0; i<32; i++) 
  64.     { t=lehmer()+lehmer(); /* lehmer is only 31 bits */
  65.       for(j=0; j<32; j++)AMem[j]+=(t>>j&1)<<i; /* put bit t{j} in bit
  66. AMem[j]{i} */
  67.     }
  68.  }
  69.  
  70.  
  71. You must call the initAMem() function once before using alea(). You can
  72. eliminate the
  73. calls to lehmer() and supply any 32-bits non-zero values, provided they
  74. are all different
  75. and contain no obvious pattern. 1,2,3,4, .. 32 are horrendous startup values.
  76.  
  77. The alea() function maintains 32 LFSRs in lockstep. Each bit of the output
  78. longword is
  79. the output bit of a LFSR keyed with a different startup value, so each bit
  80. follow a
  81. different sequence than the others.
  82.  
  83. For maximum speed, use alea() in a short loop, to ensure that the AMem[64]
  84. array stays in
  85. the on-chip cache. Moreover, avoid casting the alea() result to a float
  86. for the common practice
  87. of scaling it into the 0.0 to 1.0 range. That operation will cost far more
  88. than the whole
  89. call to alea(). Even a simple modulus "1+alea()%100" is as costly as
  90. alea() itself.
  91.  
  92.  
  93. int vec[100000];
  94. for(int i=0; i<100000; i++) vec[i]=alea()%1024 +1; /* recommended use */
  95.  
  96. On my slowpoke 25 Mhz 68040 Mac, I can make about half a million calls to
  97. alea() by second;
  98. 543261 calls/sec, last benchmark. In the same benchmark, the lehmer()
  99. function, which is typical
  100. for a linear congruential RNG, averaged 144313 calls/sec. So alea() is
  101. three to four times faster
  102. than the RNG you use. Likely even faster.
  103.  
  104.  
  105.  "JAMAIS un coup de dΘs n'abolira le HASARD, 
  106.    fut-il roulΘ dans des circonstances Θternelles." (StΘphane MallarmΘ)
  107.  
  108. Michel Saucier. michel_saucier@cmq.qc.ca
  109.